Numerical computation of solutions of the critical 
nonlinear Schrodinger equation after the singularity 



Panos Stinis 
Department of Mathematics 
University of Minnesota 
Minneapolis, MN 55455 



Abstract 

We present numerical results for the solution of the ID critical 
nonlinear Schrodinger with periodic boundary conditions and initial 
data that give rise to a finite time singularity. We construct, through 
the Mori-Zwanzig formalism, a reduced model which allows us to fol- 
low the solution after the formation of the singularity. The computed 
post-singularity solution exhibits the same characteristics as the post- 
singularity solutions constructed recently by Terence Tao. 



Introduction 

Nonlinear Schrodinger equations with power nonlinearities have been an 
intense subject of research both analytically and numerically (see e.g. [161 
[T7] and references therein). Depending on the sign and the order of the 
nonlinear term, the solutions of such equations can exhibit a varied range 
of behaviors from scattering to solitons to finite time singularities. In the 
current paper we are interested in the critical Schrodinger equation in one 
spatial dimension. It has been shown that given large enough initial data, the 
solution of this equation can exhibit finite time blow-up [111 fH] . We focus 
here on the ID case because it facilitates the numerical analysis and because 
the numerical results we present are novel even in this case. However, our 
techniques can be applied to cases with more spatial dimensions and different 
degrees of the power nonlinearity. 

We are particularly interested in investigating numerically the behavior 
of solutions to the equation after the formation of the singularity. There 
exists a large body of numerical research on the behavior of the solutions as 
they approach the blow-up instant [16]. Yet, to the best of our knowledge 
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there is no prior numerical work done on what happens to the solution after 
the singularity and whether it is possible to follow the solution after the 
singularity has formed. 

The reason that makes the computation of the solution at the singularity 
(and possibly after) difficult, is that the solution, loses its smoothness at 
the singularity. It may lose its smoothness for later too (as it happens 
with shocks) but even the loss of smoothness at one instant is enough to 
cause severe numerical problems. Loss of smoothness means that there is 
propagation of activity down to the zero length scale. Thus, the question 
is how can one represent accurately such a solution since any numerical 
calculation can only afford to resolve a finite number of length scales (called 
resolution hereafter). 

We have chosen to address this issue through dimension reduction (see 
e.g. [7] for a review) and in particular the Mori-Zwanzig formalism [3]. The 
main idea behind dimension reduction is to divide the available resolution 
into resolved and unresolved variables (length scales in our case) and con- 
struct a model for the resolved variables. The remaining computational 
capability is used to transfer activity from the resolved to the unresolved 
variables. The effect of the unresolved variables on the model equations for 
the resolved variables is to enhance these equations by terms which account 
(inevitably in an approximate manner) for the interaction between resolved 
and unresolved variables. In essence, what one looks for is a way to sim- 
plify the dynamics, since a finite resolution does not allow keeping all the 
dynamics, while at the same time retaining the most important features 
of the dynamics. As is expected [5], the main difficulty in constructing a 
reduced model is to estimate the correct rate at which activity is trans- 
ferred between the resolved and the unresolved variables. For the nonlinear 
Schrodinger equation we have chosen to use a reduced model stemming from 
the Mori-Zwanzig formalism called the i-model (see [U EJ |9] for thorough 
discussions and other applications of this model) . 

Recently, Tao p2] (see also [U EJ H2]) has constructed solutions for the 
critical nonlinear Schrodinger equation in at least four spatial dimensions 
which hold even after the formation of the singularity. The main feature of 
these solutions is that they eject a finite amount of mass instantaneously at 
the singularity instant. These solutions have at most a finite number of mass 
ejection events (depending on the magnitude of the initial data). Between 
these events the solution continues to conserve mass. For the special case 
of spherical symmetry, the amount of mass ejected is at least that of the 
ground state mass. The solutions of the i-model in ID do eject mass in a 
very narrow time interval around the singularity instant. The amount of 
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mass ejected is smaller than the ground state mass. It is not known yet 
rigorously how much mass should be ejected in the ID case (see Theorem 3 
and related comments in [12]). 

In addition to mass ejection at the singularity, the solution can lose its 
regularity for later times [TS1 Q2] • Bourgain and Wang pQ have constructed 
solutions to the ID critical Schrodinger with algebraic blow-up rate which 
retain their smoothness after the singularity. However, these solutions are 
unstable to perturbations and have not been observed numerically. There 
exist also more stable blow-up solutions (with a log-log correction to the 
blow-up rate) [141 [T2] which do lose smoothness after the singularity. The 
initial condition we have chosen gives rise to a blow-up solution according 
to the log-log scenario. This was confirmed by an independent calculation 
based on a mesh refinement algorithm developed recently by the author |15| . 
Thus, we expect a good reduced model to be able to show the ensuing loss of 
smoothness after the singularity. Even though any numerical calculation has 
finite resolution, the reduced model should behave in a way that is consistent 
with the theoretical predictions. Indeed, the solution of the i-model after 
the singularity shows an increasing roughness with increasing resolution. 
This trend suggests that the solution does indeed lose its regularity after 
the singularity. 

The paper is organized as follows. Section [1] presents some generalities 
about the ID critical Schrodinger equation. In Section [2] we give a very brief 
presentation of the i-model (see |9j for more details). Section [3] contains the 
numerical results. In Section [3. II we examine the evolution with time of the 
mass, the Hamiltonian and the li norm of the gradient of the solution. In 
Section 13.21 we examine the smoothness of the solution after the singular- 
ity. Finally, Section [J] concludes with a discussion of the results and some 
directions for future work. 



1 The critical focusing nonlinear Schrodinger equa- 
tion 

The ID critical focusing Schrodinger equation |16j is given by 

i— + Au+ \u\ 4 u = 0. (1) 
at 

The equation needs to be supplemented by an initial condition u(x, 0) = 
uq(x) and boundary conditions. We solve (jTJ in the interval [0, 2tt] with pe- 
riodic boundary conditions. If we assume that the solution remains smooth, 
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it is straightforward to show that the solution of ([I]) conserves the mass 
M{t) given by 



M(t) 



[0,2tt] 



\u(x, t)\ 2 dx 



as well as the Hamiltonian H (t) given by 

1 



H(t) 



[0,2tt] 



ol Vu(x,t)\ 2 - -\u{x,t)f 



dx, 



The use of periodic boundary conditions allows us to expand the solution in 
Fourier series 

u K ( x ,t)= Mty kx , 

fcG-FUG 

where F U G = [— ^ e have written the set of Fourier modes 
as the union of two sets in anticipation of the construction of the reduced 
model comprising only of the modes in F = [— y , y — 1], where N < K. 
The equation of motion for the Fourier mode u k becomes 



du k .,2 
at 



E 



k!...,k 5 £FUG 

where u* k denotes the complex conjugate of the Fourier mode Uk- The ODE 
system ([2]) conserves the discrete versions of the mass and the Hamiltonian 



M(t)= \ u k(t)\ 2 and H(t) = ^ 

fcGFUG keFUG 



\k 2 \u k {t)\ 2 -\\u k {t)f 

2 D 



2 The i-model 

The solution of ([T]) can blow-up in finite time depending on the magnitude 
of the initial condition uq(x). The representation of the solution of ([T]) by 
a finite number of Fourier modes breaks down at the blow-up instant since 
the solution develops activity down to the zero-scale. This presents a major 
problem for numerical calculations since we can only afford a finite number 
of Fourier modes. In other words, no matter how large a calculation we 
can afford, we are bound to run out of resolution at the blow-up instant. 
However, in some cases, it is possible to construct a model for a reduced set 
of Fourier modes (called the resolved modes) which remains well-resolved 
even after the blow-up instant. Such a model needs to be able to eject mass 
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at the correct rate from the resolved to the unresolved modes. As expected, 
the main difficulty in constructing the reduced model lies in estimating the 
correct rate of mass ejection from the resolved to the unresolved modes. 
In general, the problem of estimating the rate at which activity propagates 
from the large scales to the small scales of the solution is hard (see [7] for 
extensive discussions and examples). 

The ODE system ([2]) for the modes F L) G can be rewritten as 

= mt)), 

where u = ({itfc}), k E F U G is the vector of Fourier coefficients of u. Also, 
R(u(t)) is the vector of right hand sides (RHS) of (j2J). The vector of Fourier 
coefficients can be written as u = (u,u), where u are the resolved modes 
(those in F) and u the unresolved ones (those in G). Similarly, for the RHS 
vector we have R(u) = (R(u),R(u)). 

We need to choose a reduced model for the modes in F. We use a reduced 
model, known as the t-model, which has been shown to follow correctly the 
behavior of the solution to the inviscid Burgers equation even after the 
formation of shocks [21 [9] . It has also been used to investigate the possible 
finite time blow-up for the 3D Euler equations of fluid mechanics [9j. The 
i-model was first derived in the context of statistical irreversible mechanics 
using the Mori-Zwanzig formalism [3] and was later analyzed in [21 [9]. It 
is based on the assumption of the absence of time scale separation between 
the resolved and unresolved modes. For a mode u' k in F the model is given 
by 

U k = -ik U k +l 2^ U fci U fc2 U *8«*4 U *5 



dt 



ki...,k 5 eF 



+ t 



3 * Yl R ki{u')u'k 2 u' k3 u% A u' ks 



ki—k-2+k3—k4+k5=k 
fciGG,fc 2 ,...,fc 5 eF 



k 1 -k 2 +k 3 -k 4 +k 5 =k 
k 2 £G,k 1 ,k 3 ,...,k 5 £F 



■ (3) 



where we have suppressed the dependence of the solution on time to avoid 
clutter in the formulas. Also, the prime is used to denote the fact that the 
solution of the reduced model (|3j) can differ from the solution of the system 
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([2]). Of course, we hope that the reduced model will be able to reproduce 
the correct behavior for the resolved modes. 

Note that the t-model is closed in the resolved modes. The second term 
on the RHS of (J3]) is of the same form as the second term in ([2]), except that 
the term in ([3D is defined only for the modes in F. The third and fourth 
terms in ([3]) are not present in ([2]). They are of order nine in the Fourier 
modes and they are effecting the drain of mass out of the modes in F. The 
name of the i-model comes from the explicit dependence on time of the RHS 
of the system ([3]) . 

The i-model has two features which make it attractive as a reduced 
model for problems without time scale separation over a large range of modes 
(scales). The first feature is that it is derived directly from the system of 
ODEs and does not involve adding regularizing terms by hand. This facili- 
tates the proof of its convergence with increasing resolution to the solution of 
the PDE as long as the solution is smooth. Also, it does not involve any ad- 
justable parameters in contrast to artificially added regularizing terms. The 
second welcome feature of the i-model is that, for systems of ODEs which 
conserve the 1% norm of the solution (the mass in the nonlinear Schrodinger 
case), the I2 norm for the resolved modes is non-increasing in time [9]. In 
fact, for Mp{t) = YlkeF \ u 'k(^)\ 2 we nave 



We should note that the i-model has a similar form for the critical non- 
linear Schrodinger in more than one spatial dimensions. It can also be 
constructed for the supercritical Schrodinger equation (in one or more spa- 
tial dimensions). We have applied the t-model in these cases and we will 
present those results in a future publication. 

3 Numerical results 

In the numerical experiments we used the (purely imaginary) initial condi- 
tion 



For this initial condition we have max|ito(x)| = A at x = ir. We present 
results for the case of A = 1.80 which was found through a mesh refinement 
algorithm |15| to lead to a finite-time blow-up (more details on the behavior 
at blow-up are given later). 





keG 



uq(x,0) = iAex.p(— [x — it) 2 ). 
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Before we present the numerical results we should comment on the choice 
of the range of the resolved modes (those in F) and the unresolved modes 
(those in G). As can be seen from the RHS of the reduced model ©, for the 
critical Schrodinger in one spatial dimension, the f-model allows interactions 
of a mode with modes which have at most five times larger wavenumber. 
This means that the optimal division of an available resolution (available 
number of Fourier modes for a simulation) is to take the resolved range as 
one fifth of the available range. For example, if we can afford to calculate 
with, say K = 50 Fourier modes, we should construct a reduced model for 
N = 10 modes and use the other 40 modes as unresolved. This would mean 
that F = [-5, 4] and G = [-25, -6] U [5, 24]. In the figures, the number N of 
modes denotes the number of resolved modes. For example, iV = 256 means 
that the calculation of the t-model involves 256 * 5 = 1280 modes. The form 
of the i-model in Fourier space makes it possible to use FFT to calculate 
the nonlinear sums. If we opt for the division of modes just discussed above, 
then all the FFTs are dealiased by construction (see [9] for more details). 

The largest resolution we have used is N = 512 which means that the 
calculation involves 2560 Fourier modes. At first sight this may seem a small 
resolution for a ID calculation. However, the situation is more complicated. 
We have used the Runge-Kutta Fehlberg method of order 4-5 to integrate 
the equations of the t-model with the error control tolerance set to 10 -10 
[8]. This leads to a stepsize of about 10 -5 before the singularity. As is 
known, the error control depends on the magnitude of high order temporal 
derivatives. So, at the singularity, where the solution changes very rapidly 
in time, the stepsize will be decreased. Indeed, the stepsize falls to around 
10~ 7 shortly before and after the singularity and then plateaus to about 
10~ 6 for the remaining calculation. If one wants to use larger resolution, 
even in ID, it is advisable to parallelize the algorithm. 

3.1 Mass and Hamiltonian evolution 

Figure [1] shows the evolution of the mass for the resolved modes for dif- 
ferent resolutions, from N = 16 to N = 512. The behavior of the mass 
evolution leads us to a few observations. The first observation is that as 
we increase the resolution the mass of the resolved modes remains constant 
for a longer time. We know from theory that the mass of the solution of 
the critical Schrodinger equation is conserved for all times that the solu- 
tion remains smooth. However, the moment the solution blows up, the 
solution loses its smoothness and there is no reason why the mass should 
be conserved. In fact, Tao [18] (see also [I]) has constructed solutions for 
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Figure 1: Evolution of mass for resolved modes. The vertical line denotes 
the numerically estimated blow-up instant T = 0.13504 calculated with a 
mesh refinement algorithm. 

the critical Schrdodinger equation which at the blow-up instant eject mass 
instantaneously and then continue with a constant (but lower) mass. 

In addition, Tao showed that the number of mass ejection occurrences 
is finite. While in the current example we have found only one mass ejec- 
tion occurrence we have conducted numerical experiments with larger initial 
conditions which give rise to more than one mass ejection occurrences, but 
always a finite number of them. This is in direct contrast to the behavior of, 
say, the inviscid Burgers equation where once a shock has been established, 
the loss of mass (called energy in the relevant literature) is not instantaneous 
but persists in time and in fact it follows a power law |10j . 

The second observation is that the time of occurrence of the mass ejection 
as predicted by the t-model is in remarkable agreement with the estimated 
blow-up time based on a mesh refinement algorithm developed earlier by 
the author [15]. In other words, the t- model kicks in only when needed and 
does not eject mass unnecessarily. For our initial condition, the algorithm in 
|15j shows that the solution blows up following the log-log scenario proven 
by Perelman [14j . In fact, 

|V.(x, t )| fa ~( ln "^ t - t " ) i for t ^T-, 
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where T is the blow-up instant. For A = 1.80, the estimated blow-up time 
T = 0.13504. We monitored the mass dissipation rate of the t-model given 
by equation above and found it to be sharply peaked at T* = 0.13508, 
which also corroborates the agreement of the t-model behavior with the 
occurrence of a singularity predicted by the mesh refinement algorithm. 

The third observation is that the value of the mass after the singularity, 
as predicted by the t-model, appears to be converging as we increase the 
resolution of the reduced model. Of course, one example is not enough to 
infer the converging properties of the t-model especially since we do not 
know if it converges to the right solution after the singularity. 

A related issue is the amount of mass concentrated at the blow-up point 
and how much mass is ejected at the singularity instant. In [13J it was shown 
that in ID and for periodic boundary conditions, the solution concentrates 
at the blow-up point mass equal to the mass of the ground state on the 
whole line. The ground state equation is given by [16] 



The mass of the ground state Mq = f R \Q(x)\ 2 dx was found to be approxi- 
mately 2.7412. 

We calculated the mass concentration around the blow-up point which 
is at 7r. We find that at the instant of the singularity, a mass amount equal 
to the ground state mass is concentrated in the region [it — 0.05, it + 0.05]. 
Of course, since we can only afford a finite resolution, the area in which the 
ground state mass is concentrated is not a single point but lies in a narrow 
range around the blow-up point. After the singularity has occurred we find 
that the mass ejected by the t-model (with N = 512) is about 0.446. The 
amount of mass ejected is small compared to the ground state mass. Merle 
and Raphael [12] talk about radiative mass ejection for the log-log blow-up 
scenario, however, at this point we do not know how large this mass ejection 
should be (see Theorem 3 in |12] and the related comments) . In |18] , Tao has 
shown that under spherical symmetry and in at least four spatial dimensions 
the amount of mass ejected should be at least the mass of the ground state. 
However, it is not known yet whether these results apply to lower spatial 
dimensions too. 



AQ + |Q| 4 Q - Q = 0. 



In ID, its solution on R is given by 
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Figure 2: Evolution of gradient I2 norm for resolved modes. The vertical line 
denotes the numerically estimated blow-up instant T = 0.13504 calculated 
with a mesh refinement algorithm. 
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Figure 3: Evolution of Hamiltonian for resolved modes. The vertical line 
denotes the numerically estimated blow-up instant T = 0.13504 calculated 
with a mesh refinement algorithm. 
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Figure 4: Magnitude of the solution of the i-model with N = 512 at different 
instants before and after the singularity (estimated singularity time is T = 
0.13504). 

Figure [2] shows the evolution of the I2 norm of the gradient for the 
resolved modes, i.e. ^2 k£F k 2 \u' k \ 2 . As we have already said the L2 norm of 
the gradient of the solution should blow up at the singularity. The behavior 
of the numerically computed I2 norm of the gradient for the resolved modes 
is consistent with such a behavior and is again in remarkable agreement 
with the estimated time of the singularity. Similarly, as seen in Figure 
OH the value of the Hamiltonian has a jump of increasing magnitude with 
increasing resolution and the instant of the jump agrees again very well with 
the estimated blow-up instant. 

3.2 Smoothness of solution after the singularity 

Our purpose in this section is to provide more details which show that 
the solution computed from the i-model before, and especially after the 
singularity, is consistent with the theoretical results about the log-log blow 
up scenario in |12j and the post-singular solutions in [18]. In particular, as it 
is shown in Merle and Raphael's work (see Theorem 3 in |12j). for the log- log 
blow-up scenario, which our solution also follows, the quantization of mass 
at the singularity leads to a decomposition of the solution in a component 
that blows up in a self-similar fashion and another component, which is in 
L2 but not in H . In this case [18], Tao's post-singular solutions exist but 
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they are only in L2, i.e. they have lost regularity. Inevitably, in a numerical 
simulation we are bound to have finite resolution. Yet, we want to show that 
the trends exhibited by the i-model solution as we increase the resolution 
are consistent with the theoretical results. 

As we have seen in Figure [51 the gradient I2 norm seems to increase 
without bound (as it should) at the singularity as we increase the resolution. 
More importantly, its value after the singularity also increases with increased 
resolution. We have to note that the value of the gradient I2 norm after the 
singularity is not a constant even though it may appear to be so from the 
figure. In fact, it oscillates mildly (about 0.04%) around a well-defined mean. 

Even more abrupt is the increase with resolution of the constant value 
of the Hamiltonian after the singularity (see Figure [3]) . Since the mass value 
after the singularity appears to have converged with resolution (see Figure 
[I]), the abrupt increase of the Hamiltonian with resolution must be due 
to the increase of the gradient I2 norm. If the trends shown in Figure [2] 
continue with even higher resolution, then the gradient I2 norm should not 
only blow-up at the singularity but remain infinite after the singularity too. 
That would mean that the solution has lost its regularity. 

Figure H] shows the magnitude of the i-model solution with N = 512 at 
different instants before and after the singularity. Even though the evolution 
leading to the singularity involves the concentration of the mass at one spa- 
tial point, long after the singularity has happened and mass has been ejected 
out of the resolved modes, the solution appears to be scattering towards the 
boundary of the spatial domain. We want to see how the smoothness of the 
solution after the singularity changes as we increase the resolution. Figure [5] 
shows the form of the solution for different resolutions at time t = 0.75. We 
see that as the resolution is increased the solution becomes rougher. This 
increase in roughness with resolution is consistent with the solution losing 
its regularity after the singularity. 

To support our claim about the increase in roughness with increasing 
resolution we also look at the mass spectrum. In Figures [6] and [8] we show the 
spectrum of the t-model solution for different resolutions at time t = 0.1355 
and t = 0.75 respectively. We make two observations. First, the i-model 
solution even after the singularity (t = 0.1355), remains well resolved. For 
the highest resolution here, N = 512, the 80 highest wavenumbers have a 
total mass of O(10" 5 ) while the total mass of the solution is 0(1). Second, 
the spectrum becomes shallower as we increase the resolution. 

The previous statement can be made more precise. In Figure [7] we have 
plotted the total mass of the resolved modes with wavenumber greater or 
equal to 25 shortly after the singularity (t = 0.1355). The choice of the 
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Figure 5: Magnitude of the solution of the t- model with increasing resolution 
at time t = 0.75. 




Figure 6: Linear- log plot of the mass spectrum of the solution of the t- model 
with different resolutions at time t = 0.1355. 
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Figure 7: Total mass in wavenumbers greater or equal to 25 as a function 
of the resolution at time t = 0.1355. 

value 25 is because for the largest resolution here (N = 512), this is the 
wavenumber where the linear tail (in linear-log coordinates) of the mass 
spectrum begins. We see that as the resolution is increased there is an 
increase in the total mass in the modes with wavenumber greater or equal 
to 25. The total mass in these modes as a function of the resolution can be 
fit by a straight line with a correlation coefficient of about 0.998. Since the 
total mass in all the modes appears to converge with increasing resolution 
(see Figure [1]), the increase in mass for the wavenumbers larger or equal to 
25 must be attributed to a decrease in the mass of the wavenumbers less 
than 25. In other words, as the resolution is increased the large scales of 
the solution have a decreasing mass. Consequently, the solution becomes 
rougher. 

For t = 0.75, (see Figured]) the total mass of the 80 highest wavenumbers 
is O(10 -14 ) and the spectrum again becomes pronouncedly shallower with 
an increase in the resolution. A similar analysis as in the case shortly after 
the singularity shows that the total mass in the wavenumbers greater or 
equal to 25 also increases with increasing resolution. So, whether shortly or 
long after the singularity, as the resolution is increased the solution becomes 
rougher. 
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Figure 8: Linear-log plot of the mass spectrum of the solution of the t-model 
with different resolutions at time t = 0.75. 

4 Discussion 

We have presented numerical results from the application of a reduced 
model, called the t-model, to the ID critical Schrodinger equation for an 
initial condition that leads to a finite time singularity. The results suggest, 
in agreement with the recent theoretical results in [18] , that the singularity 
causes a finite amount of mass to be ejected instantaneously, after which, the 
solution continues to conserve the remaining mass. The time of occurrence 
of the singularity as predicted by the t-model is in very close agreement to 
the estimated singularity time given by a mesh refinement algorithm. It is 
very encouraging that the reduced model calculations are consistent with 
the blow-up instant predicted independently through a mesh refinement al- 
gorithm. Such consistency means that reduced models can help to shed light 
in open problems about the behavior of time dependent partial differential 
equations with singular (or near-singular) solutions. 

In addition to the instantaneous mass ejection, the i-model solution's 
roughness increases with increased resolution. This trend suggests that the 
solution loses regularity after the singularity. This is again in agreement 
with the theoretical results |12| . 

Finally, similar behavior for the solution to the one shown in the current 
paper has been obtained for solutions of the critical Schrodinger equation 
in dimensions two and three. As we have already seen, even in one dimen- 
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sion, as we increase the resolution the calculation becomes rather expensive 
when the solution passes through the singularity. The calculations for the 
cases with more spatial dimensions are also very expensive as we increase 
the resolution and a detailed study definitely requires parallelization of the 
algorithm. This is ongoing work and the results will be presented elsewhere. 
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